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A Landau fluid model for a collisionless electron-proton magnetized plasma, that accurately re- 
produces the dispersion relation and the Landau damping rate of all the magnetohydrodynamic 
waves, is presented. It is obtained by an accurate closure of the hydrodynamic hierarchy at the level 
of the fourth order moments, based on linear kinetic theory. It retains non-gyrotropic corrections 
to the pressure and heat flux tensors up to the second order in the ratio between the considered 
frequencies and the ion cyclotron frequency. 
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I. INTRODUCTION 

In many spatial and astrophysical plasmas, collisions are negligible, making the usual magnetohydrodynamics 
(MHD) questionable. The presence of a strong ambient magnetic field nevertheless ensures a collective behavior of 
the plasma, making a hydrodynamic approach of the large-scale dynamics possible and even advantageous, compared 
with purely kinetic descriptions provided by the Vlasov-Maxwell (VM) or the gyrokinetic equations. It is thus of 
great interest, both for the numerical simulation of broad spectrum phenomena and for an easier interpretation of the 
involved processes, to construct fluid models that extend the MHD equations to collisionless situations by including 
finite Larmor radius (FLR) corrections and Landau damping. In a fluid formalism, FLR corrections refer to the part of 
the pressure and heat flux tensors associated with the deviation from gyrotropy. They play a role when the transverse 
scales under consideration extend up to or beyond the ion Larmor radius (fluid models are always limited to parallel 
scales large compared with the ion Larmor radius). Evolving on a shorter time scale than the basic hydrodynamic 
fields, FLR corrections can generally be computed perturbatively. This expansion cannot however be pushed arbitrary 
far and any fluid analysis addressing (transverse) scales comparable to the ion Larmor radius^ can only be heuristic. 

From Vlasov equation it is easy to derive a set of exact moment equations. This fluid hierarchy is however faced 
with a closure problem. An interesting approach consists in closing this hierarchy by using relations, derived from 
linearized kinetic theory, between given moments and lower order ones. This in particular accounts for linear Landau 
damping in a fluid formalism. Such an approach initiated in Ref. [2] leads to descriptions usually referred to as 
Landau fluids. We here concentrate on a closure at the level of the fourth order moments, which provides an accurate 
description of most of the usual hydrodynamic quantities. 

An alternative method to the Landau fluids is provided by the gyrofluids 3,4 obtained by taking the moments of 
gyrokinetic equations. The same closure problem is encountered for the moment hierarchy. The gyrofluids have 
the advantage of retaining FLR corrections to all order relatively to the transverse scale within a low frequency 
asymptotics but, being written in a local reference frame, the resulting equations are more complex than those 
governing the Landau fluids, we are here concerned with. 

As an example, Landau fluid models should be most useful to analyze the dynamics of the magnetosheath that 
appears as a buffer between the earth bow shock and the magnetopause and plays an important role in decreasing 
the impact of solar activity on the earth environment. Recent analyses of data provided by the Cluster spacecraft 
mission have revealed that the magnetosheath displays a wide spectrum of low frequency modes (Alfven, slow and 
fast magnetosonic, mirror whose wavelengths extend down to the ion gyroradius and beyond. Since the plasma 
is relatively warm and collisionless, Landau damping and FLR corrections are supposed to play an important role. 
Coherent solitonic structures (magnetic holes and shocklets) are also observed, and their origin is still debated**^ 

A Landau fluid model for collisionless purely magnetohydrodynamic regimes^ was first derived from the equation for 
the distribution function of the particle guiding centers, taken to lowest order. It is thus restricted to the largest MHD 
scales where the pressure and heat flux tensors for each species can be viewed as gyrotropic and where the transverse 
velocity reduces to the electric drift. Starting directly from the VM equations, this model was then extended in order 
to include a generalized Ohm's law and to retain the leading order FLR corrections to the pressure tensorAifi This 
model enabled one to reproduce the dynamics of dispersive Alfven waves propagating along the ambient field both in 
the linear and weakly-nonlinear regimes and to recover the kinetic derivative nonlinear Schrodingcr (KDNLS) equation 
in a long-wave asymptotic expansion with, as the only difference, the replacement of the plasma response function 
by its two or four poles Pade approximants. It also accurately describes the dissipation of oblique magnetosonic 
waves>ii Non-gyrotropic contributions to the heat fluxes were introduced in Ref. [12] in order to obtain the dispersion 
relation and the Landau damping rate of oblique and kinetic Alfven waves. The approach we present here provides 
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a more systematic description of the FLR corrections up to second order, by retaining parallel and transverse heat 
flux vectors whose coupling to the non-gyrotropic pressure contributions is in particular required for an accurate 
description of the transverse magnetosonic waves^ A recent paper by Ramos^ addresses a similar issue and derives 
a complete set of nonlinear equations for fluid moments up to the heat flux vectors, leaving the closure on the fourth 
order moments unspecified. We here follow a similar path choosing in Section II to linearize the equations for the 
("slaved") non-gyrotropic contributions to the pressure and heat flux tensors, while retaining nonlinear equations for 
all the other moments. While Ramos performs a first order expansion in the regime referred to as the fast dynamics 
ordering, we here keep the second order accuracy necessary for a proper description of the oblique dynamics. By 
fitting with the kinetic theory briefly reported in Section III, we also give in Section IV an explicit closure relation, 
taking into account FLR corrections, and approximating the plasma response function with four and three poles Pade 
approximants in order to recover accurate limits for Landau damping both in the isothermal and adiabatic regimes. 
As the result of such a high order approximation, one of the fourth order moments is prescribed as the solution of a 
dynamical equation. After a discussion of the resulting model in Section V, the validation of the model at the level 
of the dispersion relation of the various MHD waves is addressed in Section VI. Section VII is the conclusion where 
further extensions to a model, aimed at including a realistic description of the mirror modes, are announced. 



II. FLUID DESCRIPTION OF EACH PARTICLE SPECIES 



A. The moment hierarchy 



Starting from the VM equations for the distribution function f r of the particles of species r with charge q r , mass m r , 
and average number density n r , one easily derives a hierarchy of fluid equations for the corresponding density p r = 
m r n r J f r d 3 v, hydrodynamic velocity u r — J vf r d 3 v/ J f r d 3 v, pressure tensor p r = m r n r f{v — u r ) ® (v — u r )f r d 3 v 
and heat flux tensor q r = m r n r j(v — u r ) ® (v — u r ) <S> {v — u r )f r d 3 v, in the usual form 

d t pr + V ■ (p r u r ) = (1) 

d t u r + u r ■ Vu r + —V ■ p r - — (E + -u r x B) = (2) 

p r m r c 

r Qr I s 

d t pr + V • (u r p r + q r ) 4- p r • Vw r H — B x p r =0, (3) 

L m r c J 

where the tensor B x p r has elements (B x p r )ij = £imiB m p r ij and where, for a square matrix a, one defines 
a tr . One has {B x p r ) tr = — p r X B. In order to distinguish between scalar and tensorial pressures, bold 
letters are used to denote tensors of rank two and higher. The equation for the heat flux tensor involves the fourth 
order moment i> = m r n r J(v — u r ) ® (v — u r ) ® {v — u r ) ® {v — u r )f r d 3 v. Since at this step we are dealing with the 
various particle species separately, we simplify the writing by hereafter dropping the r subscript. The equations 
governing the heat flux elements then read 

d t qijk + vidiqijk + din jk i - ^dipimXSmiPjk + 6 mj p ik + S m kPij) 

(4) 

We here concentrate on the ion dynamics. The corresponding equations for the electrons are obtained from the 
equations for the ions by changing the sign of the electric charge (including in the cyclotron frequency) and making 
the approximation m e /m p <C 1. 



B. Pressure tensors and heat flux vectors 



In order to isolate the gyrotropic components of the pressure tensor, it is convenient to rewrite Eq. J2J for the 
pressure tensor of each particle species in the form 

pxd - 6xp = k (5) 
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where & = j— r is the unit vector along the local magnetic field and 
LB 



k = 



f| + (V-n)p + V-q+(p-Vu)* . «,) 



In this equation, _B denotes the amplitude of the ambient field assumed to be oriented in the z-direction, and Q = 

d mC 

is the cyclotron frequency of the considered particles species with charge q and mass m. Furthermore, — = dt + u ■ V 

at 

denotes the convective derivative. 

We first note that the left-hand side of Eq. J^J can be viewed as a self-adjoint linear operator acting on p, 
whose kernel is spanned by the tensors n = I — b (g> b and r = 6 (g> b. Using the symbol : to denote double contrac- 
tion, it is convenient to define the projection a of any (3 x 3) rank two tensor a on the image of this operator as 

a = a —(a : n)n — (a : t)t, which implies tr a = and a : r = 0. In particular, the pressure tensor p = P + II is 

written as the sum of a gyrotropic pressure P = p±n + pnT (with 2p± = p : n and pu = p : t) and of a gyroviscous 
stress II = p that satisfies II : n = and II : r = 0. 

A similar decomposition is performed on the heat flux tensor by writing q = S + <r with the conditions Oi^rijk = 
and UijkTjk =0. One has 

Sijk = -^{strijk + Sfriik + S^riij + Sj-rurijk + S^Tijn ik + S^r^ny) 

+Sfrjk + S'Jnk + Slnj - - \SjruTju + Sfnjnk + SfrtkTij) , (7) 

where the parallel and transverse heat flux vectors 5" and S 1 - have components Sf = qijkTjk and 2S^- = qijunjk- In 
the special case where the tensor q is gyrotropic, only the z-components qu = S» ■ b and q± = Sj_ • b are non zero. 
Transverse components are however required, for example to describe transverse magnetosonic waves^ 

We consider in this paper perturbations that are at large scale in all space directions and in time, with an amplitude 
that is relatively small. This leads us to retain the terms involving the non-gyrotropic parts of the pressure and heat 
flux tensors at the linear level only. Such an ordering implies in particular that increasing the amplitude of the 
fluctuations requires longer length scales for preserving a given accuracy. In the following, we shall thus neglect the 
cr contribution to the heat flux tensor. One indeed easily checks from the equation satisfied by <r (see Appendix 2 of 
Ref. [15]) that cr involves either nonlinear contributions or linear contributions of second order relatively to the scale 
separation parameter, and thus turns out to be negligible in the equations for the gyroviscous stress or for the heat 
fluxes, at the order of the present analysis. 



C. Dynamics of the gyrotropic pressures 

To obtain the equations for the gyrotropic pressure components, one applies the contraction with the tensors I and 
t on both sides of Eq. © to ge t 12 ! 15 

d t p± + V ■ (up±) +p±X7 ■ u - p±b ■ Vw -b + - (tr V • q - b ■ (V • q) • b 

+ ^(tr (n • Vu) s - (n • Vu) s : t + n : -p) = (8) 

-v ~ ^ dr 

d t p\\ + V • (up\\) + 2p\\ b ■ Vu ■ b + b ■ (V • q) • b + (II • Vu) s : r - II : — = 0, (9) 

which appear as the condition for the solvability of Eq. (jSJ). Note that it is important to retain the coupling to 
the gyroviscous stress (in spite of its smallness) in order to ensure energy conservation whatever the form of the 
forthcoming closure relations^ 

Since cr does not contribute at a linear level in the pressure equations, we can neglect it and write 

6 • (V • q) • b M -2(6 • S^)W ■ b + V • S 11 - 2b ■ V6 • S 11 (10) 
-(tr(V • q) - b- (V- q) • &) « V • S ± + (b ■ S' ± )V -b + b-Vb- SK (11) 
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D. Gyroviscous stress tensor 



In order to determine the non-gyrotropic contributions to the pressure tensor of the various particle species, we 
start from Eq. (J5J for the full pressure tensor. Using Eqs. ©-© governing the gyrotropic pressures, Eq. (J3J| is 
rewritten 



n x b - b x n = k + L(n) 



where 



and 



lBo 

n \b\ 



dp 

~dt 



(V • u)P + V • q + (P • \/ u y 



? + (v ■ u )n + (n ■ vuf 



(12) 



(13) 



(14) 



The elements of k rewrite 



K 



_ ]_Bo r 

13 n \b\ 



(p\\ - P-O - ^ + d k<lkij + P±(n>ikdkUj + UjkdkUi - n i:j n kl diu k ) 



+P\\(Tikd k U J + TjkdkU, - 2T tj T kl dlU k ) 



(15) 



Furthermore in Eq. i|12fl . the element of the left-hand side with ij indices reads 

e jk iU ik bi- e ik ib k Tlij = bi(e jk iTl ik + e ik iTl k j), thus suggesting a misprint in Eq. ( 3.5) of Ref. [15]. When ne- 
glecting as previously the contribution originating from <r, the heat flux term d k q k ij reduces to 



fcmn ' %j ' run 



d k q k ij w (v • .5?) = d k S k ij - -UijUmndk S k . 
In the linear approximation, we have 

(V~S), tJ = X - \d k (si + {S x -b)b t )n ]k + d k (sj- + (S x -b)b^n lk 
+bj(b'V)S} + bi(b- V)S| - 2(6 ■ VS 11 • 6)r« - -(V ■ S 1 - -b-WS 1 - • 6)n; 



(16) 



(17) 



where the derivatives act only on the heat flux components. This yields (the superscript (0) refers to equilibrium 
quantities) 



d t Ii XX - 2W. X y + P f\d X U X - dyUy) + ^(d X S£ - dyS^) = 

1 



d t U xy + 2nU xx + P y l'(d X Uy + dyU x ) + -{dySt + d^) = 



(18) 

(19) 
(20) 
(21) 



d t u xz - nu yz + P fd x u z + P \l >) d z u x + d x si + d z s\ - (pf -pf)d$ x = o 
d t n yz + nn xz + P f ) d yUz + P \° ) d zUy + dyS^ + d z sl - ( P f -p ( ° ] )d t b y = o 

together with H xx = —H yy and H zz = 0. Defining the transverse divergence of the gyroviscous stress tensor Vj_ • II j_ 
as the vector of components ^^II^j; + dyli xy , d x TL X y + d y TL yyi O^j and introducing the unit vector z in the direction 
of the ambient field, Eqs. (|18fl and (|19l) then give 

JS>) 



V± • Tlx + ^ A ^ s± x * = ~% A ^ M x ? ~ ^Q dt ( Vi ' x 



On the other hand, defining the vector Ii z = (n xz ,H yz ,H zz — 0), Eqs. (|2U|I and (|21|) rewrite 

- nn z x z + d z s[ = -VxS x - P ^v±u, - P f ] d zUl _ + (pf - p\° } )dtb ± - d t u z 



(22) 



(23) 
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E. Dynamics of the heat flux vectors 

Equation (0} for the heat flux tensor involves the divergence of the fourth order moment r, that at this step should 
be simplified in order to conveniently close the hierarchy at the present order. We first note that instead of dealing 
with the fourth order moment r, it is convenient to isolate the deviation from the product of second order moments 
by writing 

/»;, h = 1',,1'n, + P ik Pjl + PuPjk + Pijllik + /'/.II,; + PaU-jh 

+HijP lk + H lk Pji + H a P jk + pr m . (24) 

The correction term prij k i a priori includes a contribution of the form ILjII/fc + TlikTlji + TLuTLjk that we here neglect 
since, as already mentioned, contributions from the gyroviscous stress are retained in linear terms only (except in Eqs. 
© and @ in order to ensure energy conservation). This algebraic transformation allows significant simplifications 
in the forthcoming equations. Second, we make the approximation of retaining only the gyrotropic part of the tensor 
r that is then given by 

~ Hi ii ~ 

hjH = -^-(njTki + r ik Tji + TiiTjk) + r^UijTM + n ik Tji + n a Tj k (25) 
+T i:j n k i + TikUji + Tanj k ) + —^-{UijUki + n ik nji + n a n jk ). (26) 

The scalar quantities r|||| = njiknjTki, r±\\ = -rijikfiijTki and rj_± = -njikriijUkt are related to rj|||, r\\± and r±± 
(given by similar formulas with njki replaced by ryjw) by 

nm =1111-37 ( 27 ) 
n,x = n„-^ (28) 

r±x=r±±-2^. (29) 
P 

One derives the equations for the heat flux vectors by writing — = S^ k - jk H — -^—Tik and 
H J 8 dt J dt dt J 

dSf „ dTi k dS. 



Sij k J k H -T~ n jk- The first term in the above equations is given by 



dt J dt dt 

S ,^ = 2(S ,-S.,.4 + 2S .^ (30, 
and the second terms are computed using the dynamical equation for the third order moment. One gets 

. _ Sl) .?| - S J'f - i(P„P 4 , + P (t P, + PllPlt 
+Pijllki + PikHji + P a U jk + IlijPki + U lk Pji + tt a P jk )diT jk 
-(Pjl + Uji)di(^-(Tij + 2n l0 )^j - Pjidi(^-n jk U lk ^ + -Tl lk n, ]k diTlji 

-{S x ■ X7) Ui - (V ■ u)S+- - -^diUj (s^riji + S^n m jnu + Sf-riij 

+S^ n T mi n j i + S^T m in i: j + 2S\mnjkj + fltijiSfbi - -n jk dir ijk i (31) 



and 



= 2{S ± 5") -b^ + 2S^ + -(P^Pki + P k P 3l + PuP jk 
dt dt 3 dt p 

+P«n H + Pikliji + P a U jk + TLijPui + IlikPji + ttiiPjk)diT jk 
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—(Pji + Uji)di(^—(nij + STijfj - 2Pjidi[~TjkHik) + ^U ik T jk diUji 

-(S"" • V)tti - (V • u)sf - 2diu j (^S k L T :jk nu + sfrji + sfnj - -fyTubA 
+Vle ijl sfbi - Tj k dirij k i (32) 
which do not totally identify with the result of Ref. [15]. 

F. Second order approximation of the non-gyrotropic pressures and heat fluxes 

Noting by inspection of Eqs. (|3 II) and l-il'l) that the magnitude of the transverse components of the heat flux vectors 
scales proportionally to the inverse gyrofrequency of the ions, we linearize the equations for these quantities, while 
we retain the nonlinear dynamics of the longitudinal components (see Section II G). Using 9j(ri xx i + Ti yy {) — 2diT±± 
for i = x or y and d{fi Z zl — di^ux, and introducing the temperatures Tu — mpn/p and T± = mp±/p where m is the 
mass of the considered particles, one has 

T (0) » (0) 

-^V ± • n ± - fiS| x z = -2^V ± T X) - 2V_lF ±± - d t S±. (33) 
m in 



Similarly, 



T (o) (0) (0) _ (0) 

2 -^—d z n z - 05'| x z = -^V_ L TP - 2 J —T^ 0) d z b ± - Vmn - d t S l \ . (34) 

rri m. II m 



T (0) 



Combining Eqs. (12211 and 13311 and defining the square Larmor radius r, = — ^—r gives 

miV 

1 + 4r? A1W1 • II 1 = x A , u - ^—r 2 r A, V , T 



IU VI • 11 I = — Z A Li I U I U| V U| 

4 L -V - 1 ^ 217 X 2m L - 1 X - 1 



w A ^ v ^ - ^( v - ■ n - x 1+ i A ^) (35) 



2w (0) m » (0) 2 

x v^rl 1 ' - q=- r 2 L A ±U± + -z x V ± F ±± 

mi 2 2 S2 



-SfdVi'Ei-^xS 1 ). (36) 
Similarly, combining Eqs. I|23|l and (|34|l gives 



" z 



I 1 + 2 dfe £ H IIz = n x l v± ^ +^ v ^ +p\l"d z u ± - (p^ - ppa& x + dtu, 

(0) (°) _ (°) 

-TrM— V ^ T I| 1) - 2^^lfa I 6 L + Vxfp. + ft*!) (37) 
\l z \ m 11 m 11 11 / 

7^(0) _ (0) (0) _ (0) 

(1 + 2^-0 zz ) s{ = £ x (^v,if - + + ft4 

V msi 2 / 52 V m 11 m 11 

2T (0) 

~~dp 9 '(? ±s ' +p ( l 0)v ^^ +p[| 0) ^^ - (p± -p|j 0) )ft&± - ftn a ). (38) 

Note that the operators in the l.h.s. of eqs. I|35|l - (|38(l cannot be inverted for any wavenumber, indicating the limitation 
of the fluid approach to large scales, both in the longitudinal and transverse directions. At second order in terms of 

»..,/?-&■ 

n V m n 



Thk±, these equations simplify into 



P±K„ A ... P i0) „2 A V7 .^(1) 1 A.T7 r , 1 



Vx • H X = ^-z x A,. - ^A ± VxT- - + -z x ft V, • II. (39) 
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n 2 



2d (0) m » (0) 2 1 

Sf = x V^Tj 1 ' - P -j-r 2 L A ±U1 _ + -zx V ± F ±± + -z x (40) 

^ x (vj_# + p ( l ° ) Vx^ + P f ) d z u 1 _ - {pf-pf)d$x + dji i 
(o) (°) _ (°) 

-w d * (— ^ - *t^*-Tt 0) a,b ± + vxnix) (4i) 

(0) „(°) _ J°) 

p l J \\ 7 _ .. . 



Si = £ x (^VxT« - 2^^lf c^x + Vxrp. + ft5f 

-d s (v±Si + pfv±u z + p[, 0) cW - (pf -p[, 0) )a t S ± ). (42) 



971(0) 

The last term in the r.h.s. of Eq. (|39|l can be consistently replaced by —^p^ A±dtUx, that in (|40|l by 



9r) (o) T (o) T 



V±dt( — thvV A similar substitution is made in Eqs. 1411) and (1421) . the terms involving djl z and StS'l , 

TOiZ V yn°J / 

being replaced by their leading order expressions within the linear description. 

G. Simplified nonlinear equations for the longitudinal components of the heat flux vectors 

In deriving the dynamical equations governing the longitudinal components of the heat flux vectors, we retain the 
coupling to the transverse components and to the gyroviscous tensor at the linear level only, because of the presence 
of a 1/fi factor, and the assumption that the present equations are restricted to the description of the large scales. 
We retain the other couplings that include quadratic contributions with respect to the fluctuations (weak nonlinearity 
regime). Note that the variation of b z has a magnitude that scales like the square of the perturbations. One then gets 

d t Sl +V - (Slu) + 3Sld z u z + 3 Pll (b- V)(^) -px6j.-V.Lp 

+ -^(PH ~P±)d z b z + V • (f,|||6) - 3rjn V -b- {b ± • Vj_)rin = 0. (43) 
p 

Similarly, when considering the equation governing S^~, one gets 

d t Si + V • + S^V ■ u + (6 • V) - 2 P± (b ± ■ Vj_) f — ) 

+^ (d x U xz + d y Il yz ) + V • (r,| jS) + _ f x± + Fyx) (V • b) 



p / \ p j 

Rl(Pj ~ Pi 

p v - ' / \ p 

(6x • Vx)Fxx = 0. (44) 



III. LINEAR KINETIC THEORY 

Let us assume that the equilibrium state is characterized for each particle species by a bi-Maxwellian distribution 

I vr^ ( / Tfi Tn \ 1 

function fo = -, — 7-77777 — 7777 — ,„,.. ,„ exp \ — I — tttt-Vn H tttt^x ) >. For small disturbances, the perturbation f\ of the 

(27r) J ' y T ^ ^2T" 2T 

distribution function is linearly expressed in terms of the parallel and transverse electric field components that are 

conveniently written in terms of potentials, in the form E z = —d z ^ and E± = -V±$ dt A± with B = Bqz. + V x A 

c 

and the gauge condition V • A — 0. We also denote by b z the magnetic field fluctuations along the z-direction. 

The hydrodynamic moments are easily computed in a low frequency expansion, retaining only contributions 

u k z /2T (0) fcx /2T (0) T (0) fc 2 

up to order — ~ — -\ — — <C 1, with no condition on — \/ — — . Let us also introduce b = — — J- = kZrf , 
Q il V m Q y to mil 2 

( = j- — r / — and define the functions r„(6) = e~ h IJ[b) in terms of the modified Bessel function I v (b). A standard 
\kz\V 2T | [ 0) 
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calculation leads to the following results in terms of the plasma response function i?(£) = 1 + (Z(Q, where Z{Q is 
the plasma dispersion function. 

The longitudinal and transverse temperature perturbations Tn 1 ^ and T± are given by 



T, 



(i) 



T, 



(0) 



1 - R(0 + 2( 2 R(C) 



,(0) 



T, 



(0) 



( ri (6)-r (6))^-r (6) 



T 



(0) 



(45) 



and 



T 



(i) 



T 



(0) 



(^yi?(C) - l) ( - 261^(6) + 2&r (&) - r (6)) ^ - (6r x (6) - 6r (6)) J2(C) 



(0) 



-(6ri(6)-6r (6) 



T 



(0) 



i.2 



(46) 



When restricted to the linear approximation, the elements of the heat flux tensor reduce to 
Qijk = n^m J ViVjVkficPv - UiPjjj - Ujpf^ - Ukpfj 1 ■ For the flux vectors Sf = qijkbjbk and = ^q tj k{Sjk ~bjb k ), 
one then has S 1 } = n^m J vtvjf^v - pj 0) (uj + 26 i3 u z ) and = J v^J^v - pf(2 Ui - 5 a u z ). 



It results that 



sl! = -pf, 0) ^ 



t\ 0) u 



I T (o) k 2 



(i - 3i?(0 + 2( 2 R(o) [(r (6) - rx(6)) |l + r (6) 



(0) 



(47) 



and 



± X T (°) k z 



2br Q (b) - r (6) - 2bv 1 (b))R(0-^- + -Hr„(fc) - r x (&) #(c) 



So 



T, 



(0) 



(0) _ „(0) 



J^^f b r (b)-br 1 (b) s ' 



T 



(0) 



(*-*)}■ 



(48) 



For a gyrotropic equilibrium distribution function, symmetric in the direction of the ambient field, the el- 



ements of the fourth order moment perturbation read r, 



(i) 

ijkl 



i(0), 



J ViVjVkVifid 3 v. One computes the scalar 



quantities rjjj? = r^} k bj>jb k bi — n^m f v^fid 3 v, rf^ = i^ife^w - hb^buh = ^n^m j vfvj_fxd 3 v and r^ 1 } 
I r Mfe(% _ bibj)(Si k - bkh) = \n^m J v\fid 3 v. After linearization of Eqs. l{2"T|) - l|2*9"|l one gets 

p{i 0) t[ 0) 



2C 2 (i + 2C^(0) +3(i?(C)-i) -i2^i?(C) (ri(6)-r (6))^-r (6) 



T 



(0) 



ll-L 



(0)^ 



I'll = 



(0) z 



1 - i?(C) + 2C 2 R(C)) [(2bT (b) - r (6) - 26ri(6)) ^ + &(r (6) - ri(6) 
r/ \ /T (0) 

{ (4fe 4 r 1 (6) - 46 2 r (6) - &ri(6) + 3&r (6)) [-^o)R(0 - 



T 



(0) 



(49) 
(50) 



S 



+ (26 2 r x (6) + 6r x (6) - 26 2 r (6))i 1 '(C)^| T + (26 2 r (6) - 66r 1 (6)) ($ + |-($ - vp)) } 



(51) 
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IV. A LANDAU FLUID CLOSURE 

When comparing the expression of riV with those of si or T,f^ provided by the kinetic theory, one gets 



^W2C 2 (i + 2C 2 i?(C))+3(i?(C)-i)-i2C 2 i?(C) , /zrjf 

-S! = \^-^ s Sl (52) 



TO 



2Csgn(fc z )(l-3i?(C) + 2C 2 i?(C)) 



and 



_pj 0) T | f 0) 2C 2 (l + 2C 2 i?(C))+3(i?(C)-l)-12C 2 i?(C) _pW T W 
r \\\\ ~ — - 1 - /2(C) + 2C 2 i?(C) 7p = m ^7^' (53) 

One then notices that when replacing the plasma response function R by its four pole approximant 

o fn = 4-2»y7C + (8-3^)C 2 

MU 4-6 V^C + (16-9tt)C 2 +4i0?C 3 + (6tt- 16) C 4 ' 

one has the identity 

A ^ + ! "W = ^ (54) 

. , , 32 - 9tt , -2J^ mi . , , , 

with A = and it = . I his leads to the closure relation 

3tt - 8 P 3tt - 8 



T (0) ^(1) / 9T (o) ., 

II 

which identifies with Eq. (34) of Ref. [8]. Note that this closure is here established with no assumption on the 
magnitude of the transverse wavenumbers. 

On the other hand, fin can be expressed in terms of and the parallel current j z . One has 



m i = 



m 2o?(o L s « + l r ° (6) " ri(6) J7 5 M 



where va — Bq/ y 47rp(°) is the Alfven velocity and p(°) the plasma density at equilibrium. 

When dealing with rj|j_, the approximation consisting in replacing the plasma response function R by its two pole 
Pade approximant i?2 (C) = V(l — iV^C ~ 2£ 2 ), as performed to obtain Eq. (35) of Ref. [8] is not satisfactory since 

1 — -R(C) ~i~ 2C 2 -R(C) 

it does not correctly reproduce the large ( decay of the imaginary part of the fraction 2(R(() ~~ ' ^ m ^ ar 

possible overestimate of the Landau damping by Landau fluid models are mentioned in Ref. [16]. In contrast, using 



2 - iV*C _ .„„ i - R(0 + 2C 2 i?(C) 



^ 3(C) = 2-3^C-4C 2 + 2zV^C 3 ' appr0ximati0n 2^?(0 " -2 + i0FC ' ^ ^ t0 

write the evolution equation 



l-lf-^W + ^K + 'iP L 4 -0, (57) 

Vett V7r V to / 11 to L V to p /em u 'J 

where in the large-scale limit we are here concerned with, we made the expansion bTo(b) — bl \(b) « b — kj_r L . The 
notation m p is used in situations where the proton mass is to remain unchanged when turning to the corresponding 
equation for the electron. In Fourier space, the Hilbcrt transform H z reduces to the multiplication by isgnfc z . The 
convective derivative has been reintroduced to ensure Galilean invariance. 

Finally, the reduced moment r±± turns out to be totally negligible at large scales and will thus not be retained. 
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V. COMMENTS ON THE RESULTING MODEL 

The equations derived above for the ions are easily adapted to the electrons for which they greatly simplify when 
making the approximation m e /m p <C 1. This leads to neglect the non gyrotropic components of the corresponding 
pressure tensor. Note that the transverse components of the electron heat flux vectors survive due to the contributions 
of terms involving the product m e Cl e (see Section II. F). The system is to be supplemented by Faraday equation and 
Ampere's law where the displacement current is neglected. In this two-fluid formulation, energy is conserved, as 
discussed by Ramos4^ It might nevertheless be advantageous to filter out the scales associated with electrostatic 
waves by prescribing quasi-neutrality, replacing the electron momentum equation by a generalized Ohm's law, and 
turning to a one-fluid description. Numerical simulations of a monofluid model obtained from a simplified version of 
the present model have shown that energy is in practice very well conserved^ 

When compared with the previous modeU^ designed to reproduce the oblique Alfven wave dynamics, the present 
approach proves to be more systematic and, as discussed below, allows one to accurately simulate all dispersive MHD 
waves, including oblique and transverse magnetosonic waves (see Section VI). The previous model has on the other 
hand the advantage of including a nonlinear description of the gyroviscous tensor. It is of interest to see how, when 
linearized and restricted to the case of the Alfven wave scaling (also neglecting the gyroviscous tensor contribution), 
the equations governing the gyrotropic heat fluxes in the present model compare with those of the previous one. It 
turns out that Ref. [12] unfortunately includes a few algebraic errors originating from a sign error leading to an 
incorrect factor 3 in Eq. (C.8), a missing multiplicative factor m p /m r in the r.h.s. of Eqs. (C.9) and (C.10) and a 
missing minus sign in front of the first occurrence of Q, p /Q, r in Eq. (C.12). This in particular affects the equations 
for the gyrotropic heat fluxes where the contribution v\ e in the r.h.s. of Eq. (56) should be suppressed, the square 
bracket in Eq. (66) replaced by [v^sgng,. — v 2 A {5 rp — 1) — v 2 h r S rp ]/v A and the factor 3 in the last term in the r.h.s. 
of Eq. (67) also discarded. After correcting these errors and taking into account that pressure and heat flux tensors 
were computed using barycentric velocities, one easily checks that the parallel heat flux equation is exactly recovered 
and that the equations for the perpendicular heat flux of both models identify in the isothermal limit where the time 
derivatives are negligible. This limitation originates from the insufficient order of the Pade approximant used in the 
previous model. 

VI. MHD WAVE DYNAMICS 

When restricted to a one or quasi one-dimensional dynamics along the ambient field, only the longitudinal compo- 
nents of the parallel and transverse heat flux vectors (that correspond to the gyrotropic contributions to the heat flux 
tensor) arise in the equations of motion. A long-wave reductive perturbative expansion performed on the resulting 
Landau-fluid model reproduces the kinetic derivative nonlinear Schrodinger equation derived from the VM equations 
for Alfven waves with a typical length scale large compared with the ion Larmor radius^i up to the replacement of the 
plasma response function by appropriate Pade approximants. As a consequence, the modulational type instabilities 
(including filamentatior. 18 ) of Alfven waves and their weakly nonlinear developments are correctly reproduced. 9 Nu- 
merical simulations of such regimes are presented in Ref. [10] where a study of the decay instability is also presented 
and validated by comparison with hybrid simulationsii 9 - 

As stressed in Ref. [13], the correct determination of the dispersion relation for transversally propagating magne- 
tosonic waves requires a detailed description of non-gyrotropic contributions to the pressure and heat flux tensors. 
When restricted to a purely transverse dynamics, the present model reduces to the fluid model used in Ref. [13] that 
exactly reproduces the large-scale kinetic theory (note that a factor 3/2 is missing in front of the z-term in e xy given 
in Eq. (2.8) of the latter reference). 

The present model easily reproduces the dispersion relation for kinetic Alfven waves (KAW) for which the crucial 
ingredient is the contribution to the transverse velocity originating from the time derivative of the leading order 
gyroviscous stress [last term in Eq. 1|39|) ] Whereas these KAWs are also captured by a low frequency expansion 

of the kinetic equationS)22i2^ this is not the case for oblique Alfven waves. The reason is that an expansion at order 
uj/fl includes contributions of order k\r\ when k z /k± scales like k±r^ as for KAWs, but only includes terms of order 
k±r^ for finite angles of propagation. The same limitation holds for the gyrokinetic formalism. The present fluid 
formalism however allows one to obtain the correct linear dynamics for oblique Alfven waves, as was shown in Ref. 
[21], using a Landau fluid model actually contained in the present one. 



11 



VII. CONCLUDING REMARKS 



We have constructed a Landau fluid model that reproduces all large-scale dispersive MHD waves in a warm colli- 
sionless plasma. This model may be most useful not only for numerical simulations involving a broad range of scales, 
but also for analytic purposes, such as the computation of secondary instabilities. An example is provided by the 
filamentation instability of parallel propagating Alfven waves. This mechanism may be relevant in the understanding 
of the evolution of Alfven waves in magnetospheric plasmas that often display very filamentary structures^ The 
present model allows one to account for linear Landau damping, dominant FLR corrections as well as drift velocities, 
that play an important role in these plasmas whose equilibrium state often involves a large scale longitudinal current. 
The importance of nonlinear kinetic effects such as particle trapping that are here neglected should be estimated by 
comparison with fully kinetic simulations. 

In a sufficiently anisotropic plasma, the mirror instability can develop, whose threshold is accurately reproduced by 
the present fluid model. 8,10 A difficulty nevertheless originates in that, for large-scale mirror modes, the growth rate of 
perturbations propagating in the most unstable direction scales like the transverse wave number of the perturbation, 
which makes the smallest scales retained in a large-scale simulation to be the most unstable. The instability actually 
reaches a maximal rate at a scale comparable to the ion Larmor radius and is arrested at smaller scales, under the effect 
of FLR corrections^ Small transverse scales are thus to be retained. A promising approach consists in expressing, at 
the level of the linear kinetic theory, non-gyrotropic contributions in a closed form suitable for being incorporated into 
fluid equations. Explicit reference to the plasma response function should in particular be eliminated. A model that 
reproduces the arrest of the mirror instability and that is simple enough to permit accurate numerical simulations 
will be presented in a forthcoming paper 
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